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TECHNICAL  SUMMARY 

Let  A  be  the  measure  of  the  lack  of  resistance  to  sh^ck  for  a  component, 
the  life  of  which,  say  X^,  will  be  exponential.  If  the  variability  of  manu¬ 
facture  determines  the  frequency  of  the  different  A-values,  described  by  a 
r.v.,  Y  G,  then  T  =  Xy  is  the  life  length  of  a  component  selected  at  random. 

It  will  have  a  survival  distribution 

oo 

R(t)  =  EVP[XV  >  t|Y]  =  /  e"XtdG(A)  for  t  >  0. 

Y  r  0 

This  is  a  mixed-exponential  distribution  with  hazard  Q  =  -InR.  If  G  is  the 

r(a,6)  distribution  we  find  the  hazard  to  be  Q(t)  =  aln(l  +  tg),  which  is  a 
Pareto  type  II  law,  known  from  applications  in  economics. 

In  general  if  q  =  Q1  is  any  mixed-exponential  hazard  rate  then  for  t  >  0 

00 

I*e"UdG(A) 

q(t)  -  -Q= - ’ 

[  e~XtdG( A) 

J0 

and  all  such  hazard  rates  are  known  to  be  decreasing. 

Let  2  be  a  given  class  of  decreasing  hazard  rate  functions  with  the 
following  properties:  q  e  2  is  twice  differentiable,  standardized,  i.e., 
q(0)  =  1,  and  it  and  the  induced  functions^  and  5,  where 


<p(x)  =  xq(x)  and  ?(x)  =  1  +  xq’(x)/q(x) 

satisfy 

1  ip  is  increasing, 

2°  q  is  log-convex, 

3°  £  has  a  limit  at  «  and  0  <_  c(=°)  £  1, 

alternatively,  sometimes  the  stronger  condition  is  assumed,  viz., 

3'  ?  is  decreasing. 

The  unknown  parameters  of  the  life  distribution  under  study  are  introduced 
in  a  manner  consistent  with  the  gamma  mixed  exponential;  viz., 

R(t)  =  e'a^tS^ 

where  a  >s  the  shape  parameter  and  S  the  scale  parameter,  both  positive. 
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Each  element  of  this  class  D  will  generate  a  two  parameter  family, 
which  is  a  subset  of  the  DFR  distributions. 

Are  there  any  distributions  with  hazard  rates  which  satisfy  these 
conditions?  Yes,  the  exponential  and  the  gamma-mi xed-expcnential  do.  Are 
there  other  DFR  distributions  which  are  not  mixed-exponential  which  do? 

Yes,  the  failure  rate  of  the  r(v,l)  distribution  does,  when  0  <  v  <  1 
although  it  is  not  easy  to  show. 

Are  there  any  closure  properties  to  this  class?  Yes,  if  q  is  a 
decreasing  standardized  hazard  rate  satisfying  1°,  2°,  3°  (or  3')  then 
qY  for  any  0  <  y  <_  1  (called  an  Afanas'ev  generalization)  does.  Is 
2  closed  under  mixtures  of  distributions  with  hazard  rates  in  2  ? 

Not  always.  If  G  is  a  distribution  with  support  at  only  two  points,  the 
G-mixed  exponential  will,  if  the  probability  of  the  higher  failure  rate 
is  not  too  large  relative  to  the  lower. 

We  now  ask  what  are  sufficient  conditions  on  the  mixing  distribution  so 
that  assumptions  1°,  2°,  and  3°  will  be  satisfied. 

Theorem  1:  If  the  mixing  distribution  G  is  such  that  the  function 

K(x,y)  =  G(y/x)  for  x,y  >  0 

is  totally  positive  of  order  2,  then  1°  is  satisfied,  i.e.,  ip - 
is  increasing.  b 

We  now  give  some  conditions  that  a  mixed  hazard  rate  be  log-convex  on 
(0,°°) . 

Theorem  2:  If  for  a  distribution  G  on  (0,°°)  we  define  the  convolution 

G^  =  G  *  G  and  the  related  function 
s  s 

«,(.)•  I.  I  (t  -  u)2  dG(t)  dG(u)  for  s  >  0 
£  u=0  t=0 

then  the  failure  rate  of  the  G-mixed  exponential  distribution  is 
log-convex  if  and  only  if  for  every  exponential  ramdom  variable  Z 
we  have 

E  2Aq(Z)  >  EZAg'(Z), 

where 

y 

AQ(y)  =  J  sG2(y  -  s  )  dG^(s). 

0 


(1) 
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Corollary  1:  If  the  mixing  distribution  G  is  such  that  the  induced  function 
Ag,  defined  in  (1)  is  concave  increasing  then 

2Ag(x)  >_  xAg'(x)  for  all  x  ^  0 

and  failure  rate  qG  is  log-convex  and  2°  is  satisfied. 

Corollary  2:  If  the  mixing  distribution  G  has  a  density  g  *  G*  and  the  ratio 

x  g^(x)/g2(x)  for  all  x  >  0 

is  monotone  decreasing,  where  g2  and  g^  are  the  corresponding 

derivatives  of  the  functions  defined  in  theorem  2,  then  the  failure 
rate  qG  is  log-convex. 

Remark:  If  G  is  a  r(a,l)  distribution  then  one  verifies  directly  that 


which  decreases. 


x  q^(x)  _  2a  +  1 
gJCx)  “  x 


for  x  >  0 


We  do  not  consider  as  reasonable  (for  our  application)  mixing  distributions 

_x-2 

which  are  of  infinite  order,  such  as  G(x)  =  e  for  x  >  0,  at  the  origin.  We 
admit  for  consideration  only  distributions  which  are  either 

(i)  discrete  in  some  neighborhood  of  zero  or 


(ii)  of  finite  order  at  the  origin,  viz.,  there  exists 
<  >  0  such  that  x~K  G(x)  •+  c  >  0  as  x  -*■  0. 


Theorem  3:  For  any  G-mixed  exponential  distribution  the  induced  function 

CG  has  the  properties  that  Cg(0)  *  1»  is  initially  decreasing,  bounded 

above  by  unity,  and  the  limit  exists  with 

CG(«)  =  1,  if  G  is  of  type  (i)  and 
£q(°°)  =  0,  if  G  is  of  type  (ii),  moreover 
>_  0  if  and  only  if  is  increasing. 

We  now  consider  conditions  on  the  mixing  distribution  which  will  insure 
that  1  -  c  is  either  monotone  increasing  or  unimodal.  Such  behavior  can 
often  be  easily  checked  in  specific  instances,  but  we  have  a  sufficient 
condition  in 


Theorem  4:  If  the  associated  function  Ag,  as  determined  from  the  mixing 
distribution  G  in  (1),  is  such  that  the  kernel 

K(x,y)  =  AQ(y/x)  for  x,y  >  0 


0* 
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is  totally  positive  of  order  2,  then  the  function  1  - 

induced  by  the  fail ure  rate  of  the  G-mixed  exponential 
distribution  will  have  at  most  one  mode  in  [0,-»]. 

It  is  assumed  that  we  are  given  a  sample  vector  _t  =  (t^,  ...»  t^,  . ..,  t  ) 
where  t^,  . ..,  t^  are  ordered  observations  of  times  of  failure  while 
tk+1,  ....  tn  are  the  ordered  observed  alive-times  (censored  lives). 

The  two  empirical  distributions  of  the  times  of  failure  and  of  the 
sample  will  be  denoted  by 

Fk(y)  =  £  of  ti  £  y  for  i  =1,  ....  k} 

and  similarly  for  F  .  We  make  the  notational  convention,  to  be  used  sub- 

subsequently  for  any  function  g,  that  its  transform  by  an  empirical 
distribution,  is 

00  00 
?(x)  =  /  g(xt)dFn(t),  g(x)  =  J  g(xt)dFk(t). 

0  0 

Thus  the  likelihood  can  be  written 

L(o  0|t)  =  In a  +  £n6  +  £nq(3)  -  y-  Q(B) . 

Parenthetically,  for  given  6  >  0  the  likelihood  L  ( - 1 B »  t.)  is  concave  on 
(0,®)  and  the  MLE  of  a  exists  uniquely  and  can  always  be  obtained  from  the 
equation  L'(a|8*_t)  -  0.  This  tranforms  all  observations  to  an  exponential 
with  unknown  failure  rate  a  and  so  the  MLE  is  given  by 

S.  =  k/nQ(S) .  (2) 

We  examine  the  case  for  unknown  3  in 

Theorem  5:  When  the  shape  parameter  a  is  presumed  known  and  a  sample  vector  _t 
and  q  e  2  are  given,  there  exists  an  MLE  of  B>  denoted  by  B»  and  defined 
implicitly  as  the  smallest  positive  root  of  the  equation 

c(x)  -  Y  *(x)  =  0  (3) 

only  when 

inf  qf  X(y)  <  y-  <  1. 
y  >0  K 

A  simpler  situation  exists  in  the  following  case. 

Theorem  6:  If  q  e  2  is  an  Afanas'ev  generalization  of  a  G-mixed  exponential, 
when  G  is  of  type  (ii)  and  0  <  y  <  1,  with  a  known  then  for  a  given  sample 
the  MLE  of  8  exists  if  and  only  if 

1  "  Y  <  X  <  1- 


♦ 
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Corollary  3:  Under  the  hypothesis  of  theorem  6,  if  is  monotone  then  so  is  5, 

the  corresponding  function  for  the  Afanas’ev  general  ization,  and  the  MLE  of  8 
exists. uniquely,  i.e.,  there  is  at  most  only  one  solution  to  the  equation  (3). 

We  now  turn  to  the  estimation  problem  when  both  a  and  8  are  unknown. 


Theorem  7:  For  a  given  sample  _t,  with  q  e  2  specified  and  a, 8  both  unknown, 
a  MLE  of  8,  say  8,  exists  as  the  smallest  positive  root  of  the  equation 

£(x)  -  <$(x)  =  0  for  x  >  0 

where  <j>  =  ijj/Q’,  if  and  only  if  the  sample  satisfies  the  inequality 

2tt  <7.  (4) 

A 

When  8  has  been  determined,  the  MLE  of  a,  say  is  then  given  by  an 
analogue  of  equation  (2),  namely, 

a  =  k/n  Q(8)  . 

Our  computational  procedure  based  on  the  censored  sample  t_  and  the 
assumption  q  z  2  is  as  follows: 


Algorithm: 


~  _  2 

(a)  Compute  the  sample  moments  t,  t,  t  . 

(b)  If  inequality  (4)  is  not  satisfied  the  observations  are 
from  an  exponential  distribution  with  failure  rate  A.  and 
then  estimate  it  by 

X  =  . 

nt 

(c)  If  inequality  (4)  is  satisfied  then  use  the  sample  functions 
explicitly  given  by 


<f(x 


)  =  x  E  t1q(xti)/  E  Q(xt. )  c(x)  =  1  +  £  E  [t.q*  (xc.  )/q(xt. )  ] . 


We  guess  8g,  then  iterate  using  the  inductive  step; 

given  8^  compute  $(8^)  and  calculate  8^^  such  that 

c(ei+1)  =  ♦(Bi). 

We  find  8  =  lim  8^ *  and  compute 

d  =  k/  Z  Q( t. 8) •  ' 

1  1 
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TJie  nature  of  the  intersection  guarantees  that  within  the  region  when 
$  and  z  both  decrease  the  iteration  will  rapidly  converge,  with  a  reasonable 
first  guess.  When  the  functions  q  and  Q  are  simple  a  small  programmable 
electronic  calculator,  such  as  the  HP-67,  can  be  used  to  obtain  these  estimates. 

We  now  present  some  data  sets  from  two  different  lots  of  flight  control 
electronic  packages.  Each  package  has  recorded,  in  minutes,  either  a  failure 
time  or  an  alive  time,  the  latter  is  denoted  by  an  affix  +. 

First  Data  Set 

1,  8,  10,  59+ ,  72+,  76+ ,  113+,  117+,  124+,  145+,  149+,  153+,  182+,  320+. 

Second  Data  Set 

37,  53,  60+ ,  64+,  66+,  70+,  72+,  96+,  123+. 

One  checks  that  both  data  sets  satisfy  condition  (4)  so  that  both  para¬ 
meters  can  be  estimated  in  a  gamma-mixed-exponential  model.  Then  using  the 
estimation  techniques  derived  previously  we  obtain: 

Date  Set  1  Data  Set  2 

a  =  .0453,6  =  1.03  a  =  .420,  6  =  .01 

A  statistical  test  to  determine  whether  the  data  require  a  constant  or 
decreasing  failure  rate  was  run  on  the  data  from  sets  1  and  2.  For  data  set 
1  we  reject  constant  failure  rate  (in  favor  of  decreasing  failure  rate)  at 
the  .10  level.  For  data  set  2  we  cannot  reject  the  constant  failure  rate  assumption 
at  that  level.  Inthis.case,  however,  the  constant  failure  rate  estimates  for 
reliability  and  the  mixed  exponential  estimates  for  reliability  are  close. 
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Q  Abstract 

-v.. 

In  this  paper  some  of  the  problems  of  parametric  estimation  for 
two  important  types  of  decreasing  failure  rate  (DFR)  distributions  are 
discussed.  The  first  type  of  distribution  represents  the  life  of 
mixed  exponential  populations;  the  second  type  of  distribution 
represents  the  life  of  mechanisms  which  "work  harden"  as  they  age, 
i.e.,  old  components  are  better  than  new. 

A  representation  for  the  second  type  of  distribution  is  given  so 
that,  when  satisfied,  certain  functional  properties  of  the  failure  rate 
of  mixed  exponential  distributions  are  guaranteed.  Conditions  are 
then  imposed  on  the  mixing  distributions  that  will  insure  that  the 
failure  rate  of  each  type  of  distribution  satisfy  these  conditions. 

Shape  and  scale  parameters  for  any  standardized  DFR  family  of  this 
type  are  to  be  estimated  from  the  type  of  data  which  is  available  in 
practice,  namely,  severly  censored  samples  with  only  a  few  failure 
observed,  all  of  which  occur  early. 

Conditions  are  obtained  that  maximum  (rather  than  minimum)  likelihood 
estimates  (NILE’S)  exist;  the  condition  is  phrased  in  terms  of  censored 
samples,  namely  if  t^,  .  .  .  ,  t^  are  observed  failure  times,  while 
tjc+i>  .  .  .  ,  t  for  1  <_  k  £  n,  are  censored  life  observations  from 
this  class  of  DFR  distributions  then  the  NILE'S  of  shape  and  scale 
parameters  exist  if  the  following  inequality  is  satisfied: 

k  n  n 

j=l  i=l  ■-  i=l 
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"Practical  methods  for  the  computation  of  the  MLE's  are  given. 
Actual  data  obtained  from  testing  of  integrated  circuit  electronic 
packages  illustrate  the  applicability  and  utility  of  the  techniques  and 
the  results' described. 
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1 .  Introduction 

Because  of  the  demand  for  high  reliability  in  electronic 
manufacturing,  increasingly, integrated  circuit  modules  are  being  used. 

As  a  consequence,  modular  testing  now  frequently  encounters  two 
conditions:  the  first  is,  the  life  distribution  seems  to  exhibit  a 
decreasing  failure  rate  and  the  second  is,  samples  are  virtually  always 
censored. 

Of  course  these  events  are  not  unrelated,  the  first  implies  the 
second.  The  expense  of  testing  will,  when  coupled  with  a  decreasing 
percentage  of  failures  over  time  in  the  surviving  population,  nearly 
always  result  in  life  length  observations  being  censored  in  practice. 

Several  censored  life  data,  with  a  paucity  of  failures,  does  not 
allow  the  usual  statistical  methods  such  as  employing  the  Kaplan -Meier 
(1958)  estimate  of  the  survival  function  and  plotting  the  negative  of 
its  logarithms  to  see  if  its  concave  or  convex.  In  many  engineering 
applications  the  monotone  behavior  of  the  hazard  rate  function  can  be 
deduced  frcm  the  physical-chemical  nature  of  failure.  In  such  cases, 
only  parameters  need  be  estimated.  There  is  a  lack  of  estimation  pro¬ 
cedures  which  can  utilize  highly  censored  data  and  which  avoid  the 
potential  for  bias  inherent  in  Bayesian  priors. 

In  this  paper  the  maximum  likelihood  estimates  (MLE's)  are  obtained 
for  both  the  shape  and  scale  parameters  of  a  class  of  two  parameter 
families  of  decreasing  failure  rate  (DFR)  distributions  and  conditions 
are  given  for  their  existence  along  with  practical  methods  of  com¬ 
putation.  These  estimates  are  derived  for  censored  data,  which  contains 
only  a  few  failure  observations,  and  a  fortiori  for  complete  samples. 
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The  conditions  for  the  existence  of  the  MLE's  apply  to  the  case 
when  the  mean  and  variance  of  the  DFR  distributions  do  not  exist.  It 
is  thought  important  to  have  estimation  procedures  for  such  distributions. 

2.  DECREASING  FAILURE  PLATE  MODELS 

We  now  discuss  the  physical  processes  which  determine  the  length 
of  life  ’order  consideration.  Firstly  let  us  suppose  the  quality  of 
construction  of  a  component  determines  the  level  of  resistance  to  stress 
which  it  can  tolerate.  Secondly,  suppose  the  service  environment  provides 
shocks  of  varying  magnitude  to  the  component  and  failure  takes  place 
when,  for  the  first  time,  the  stress  from  an  environment  ally  induced 
shock  exceeds  the  strength  of  the  component . 

If  the  time  between  shocks  exceeding  any  specified  magnitude  is 
exponentially  distributed,  with  a  mean  depending  upon  that  magnitude, 
then  the  life  length  of  each  component  ’will  be  exponentially  distributed 
with  a  failure  rate  which  is  determined  by  the  quality  of  assembly.  It 
follows  that  each  component  in  service  will  have  a  constant  failure  rate 
but  that  the  variation  in  manufacturing  and  inspection  procedures  will 
cause  the  population  to  exhibit  a  decreasing  failure  rate.  Instances  of 
such  natural  mixing  of  exponential  distributions  were  first  discussed  by 
Proschan  (1963).  Subsequently  other  distributions  with  decreasing  failure 
rates  of  practical  interest  were  discussed  by  Cozzolino  (1968).  But 
parametric  families  of  mized  exponential  and  DFR  distributions  have 
received  little  attention  compared  with  their  IPR  counterparts. 

Alternatively,  let  us  consider  structures  subjected  to  dynamic  stresses 
of  such  a  nature  that  the  first  stresses  to  which  they  are  subjected,  if  not 
severe  enough  to  initiate  failure,  only  caused  localized  yielding  and 
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deformation,  thus  effecting  local  stress  relief  and  reinforcement.  Such 
behavior  increases  their  ability  to  withstand  future  stresses  and  could  be 
thought  of  as  "work  hardening".  In  such  cases  older  structures  in  service 
actually  have  greater  resistance  to  fatal  shocks  than  younger  ones;  i.e., 
each  component  in  service  has  a  decreasing  failure  rate.  (An -analogous 
behavior  exists  for  increased  strength  or  immunity  in  biological  systems.) 

This  may  be  thought  of  as  "the  older  the  better"  (at  least  for  certain  periods 
and  purposes).  Usually  the  failure  mechanism,  and  its  interpretation  in 
these  DFR  cases,  is  quite  different  from  that  of  a  mixed  exponential  model. 


Let  \  be  the  measure  of  the  lack  of  resistance  to  shock  for  a 
component,  the  life  of  which,  say  X^,  ’.<411  be  exponential.  If  the  vari¬ 
ability  of  manufacture  determines  the  frequency  of  the  different  values , 

which  we  describe  by  a  r.v.,  say  A,  with  distribution  G,  then  T=XA  is  the 
life  length  of  a  component  selected  at  random  from  those  manufactured;  it 
will  have  a  survival  distribution  R  which  can  be  written  as  the  conditional 
expectation 

00 

R(t)  =  EaP[Xa  >  t | A]  =  J  e'XtdG(X). 

*0 


This  is  a  mixed-exponential  (survival)  distribution.  It  is  also  a  LaPlace 
transform. 

One  important  mixed-exponential  model  is  based  on  the  gamma  density: 


If  for  some  a,  8  >  0  we  take  the  density 

a-i  -X/3 

G’(X)  =  - - - -  for  X  >  0 

T(a)  0a 

then  the  hazard  function  Q  =  -  £n  R  is 

Q(t)  =  aln(l  +  tS)  for  t  >  0 
.and  the  density  f  *  Q'e‘Q  is  given  by 


f(t)  = 


a{3 


(1  +  t3) 


a+1 


(2.0) 


'  for  t  >  0. 
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This  two-parameter  family  seems  to  have  been  introduced  in  this  country 
by  Lomax  (195*0  who  regarded  it  as  a  generalization  of  the  Pareto  distribution. 
It  was  called  the  Pareto  Type  II,  and  it  has  proved  to  be  useful  in  business 
analysis . 

The  estimation  problems  for  this  distribution  have  been  treated  in  a 
series  of  papers  by  Kulldorff  and  Varmman  (1973)  and  Vanrman  (197 6)  but 
without  the  physical  interpretation  made  here. 


Other  investigators  have  made  use  of  the  gamma  density,  as  we  have,  to 
model  variability  among  manufactured  components,  each  one  of  which  has  a 
constant  failure  rate  while  in  service,  see  the  papers  of. Harris  and  - 
Singpurwalla  ( 1963 ) ,  ( 1969 ) . 

In  certain  instances,  such  as  in  medical  studies  of  the  etiology 
of  a  disease,  when  the  hazard  rate  is  assumed  to  be  either  increasing  or _ 
decreasing  a  generalization  of  the  Pareto  law,  exhibiting  both  types 
of  behavior  depending  upon  its  parameters,  has  been  applied  by  David 
and  Feldstein  (1979) .  They  obtained  the  maximum  likelihood  estimates 
implicitly  for  two  of  three  parameters,  neither  one  being  a  scale  para¬ 
meter,  in  the  case  of  progressively  censored  data. 

The  question  we  address  is  how  should  one  proceed  when  from 
physical  circumstances  the  hazard  rate  is  known  to  be  mixed  but  by  an 
unknown  distribution? 


We  study  the  consequences  of  a  general  mixing  distribution.  Let  q  be 
any  mixed  exponential  hazard  rate,  then  for  any  t  >  0 


q(t)  = 


O 

l 


Ae'XtdG(A) 


/  6 

JQ 


-At 


(2.1) 


dG(A) 


where,  without  loss  of  generality,  the  scale  of  X  has  been  chosen  so  that 

OO 

q  is  itandaAdiz&d,  i.e.,  q(o)  =  J AdG(A)  _=  l. 


The  unknown  parameters  of  the  life  distribution  under  study  are  intro¬ 
duced  in  a  manner  consistent  with  the  gamma-mixed  exponential;  viz.,  f 

R(t)=e'^tt6)  C2.Z)  : 

where  a  is  the  shape  parameter  and  3  the  scale  parameter,  both  positive. 

We  now  introduce  a  model  for  work-hardened  DFR  life  distributions,  called  - 
an  ge.MAalUza.tion,  by  postulating  a  hazard  function  (for  unknown 

a,  3  >  0  but  known  y  with  0  <  y  <  1)  of  the  form 

tS  1 

a 

where  q  is  defined  as  in  (2.1)  as  a  mixed- exponential. 

It  is  clear  that  for  y  <  1  we  have  a  decreasing  failure  rate  distributioi 
which  is  not  mixed  exponential.  In  practice  the  value  of  y  is  often  determined 
by  the  material  properties  of  the  component  in  service  and  is  not  estimated  by 
statistical  techniques,  see  Weibull  (1961). 

Let  us  find  the  standardized  hazard  rate,  say  qA,  for  the  Afanas'ev 
generalization  of  the  gamma  mixture.  We  have 


/ 


[q  (x)]Ydx 


qA(t)  =  (1  +  t)_Y  for  t  >  0  (2.3) 

from  which  we  find  the  corresponding  hazard  function  to  be 

[(1  +  t)1_Y  -  1]/(1  -  y)‘  for  0  <  y  <  1 

(2.4) 

ln(l  +  t)  fox^  =  1 

For  y  =  1  this  is  equivalent  with  a  Weibull  distribution  with  a  location 
parameter  of  minus  one  and  a  shape  parameter  of  less  than  unity.  In  the 
Soviet  Union  this  model  was  introduced  by  Afanas'ev  (1940)  as  a  distribution 
for  fatigue  life  in  metals. 


We  now  examine  the  distribution  which  results  from  the  mixing  of  two 
exponential  distributions .  It  is  presumed  that  owing  to  occasional  laxity 


in  quality  control  there  is  a  low  probability  p  <  1/2  of  passing  a  component 
containing  a  defect  which  can  cause  a  high  failure  rate  X  >  1.  But  there  is 
a  high  probability  of  passing  a  component  having  the  nominal  (low)  failure 
rate,  which  without  loss  of  generality  we  take  to  be  unity. 

The  reliability  of  this  Bernoulli-mixed  exponential  population  is 

R(x)  =  pe-^x  +  (l-p)e-x  for  x  >  0. 

Whence,  we  find  the  standardized  hazard  to  be,  for  x  >  0 

Qg(x)  =  -£np  +  ~  -  .  £n(l  +  revx) 

One  verifies  easily  that  qB(0)  =  (pX  +  q)/u  =  1  with 

r  =  f  y  =  —A  >  0,  y  =  p(X  +  r)  . 

We  shall  write,  respectively, 

T~JA(a,B;r)  or  T  ~  JB(a,p:p,X) 

whenever  the  standardized  hazard  QA  ,  is  defined  in  (2.4)  or  Qg  as  defined 
in  (2.5).  Also  without  further  mention  we  shall  use  the  same  subscript  to 
denote  other  functions  associated  with  these  cases. 

It  is  possible  to  introduce  an  alternative  parameterization  to  (2.2), 
namely 

-  f  Q  (tp) 

R(t)  =  e  p  for  t  >  0;  a,B  >  0,  (2.7) 

where  Q  is  a  known  cumulative  hazard.  The  advantage  of  this  form  is 
that  with  standardized  q  we  see  from  L'Hospital's  rule  that  B  •*  0  implies 
R(t)  -*■  e  at.  Thus  the  limiting  case  is  an  exponential  distribution.  In 
either  formulation  one  sees  that  B  measures  the  rate  of  departure  from  a 
constant  failure  rate  in  terms  of  Q  and  hence  determines  the  rate  of 
decrease  in  the  failure  rate  with  use.  The  disadvantage  of  (2.7)  is  that 


(2.5) 

(2.6) 
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P  is  no  longer  a  scale  parameter  and  the  likelihood  equations  are  somewhat 
more  complicated.  However,  the  transformation  (a,p)  -  (a/0,0)  is  a  1-1 
mapping  of  the  positive  quadrant  into  itself  so  that  any  maximum  likelihood 
estimates  for  one  parameterization  could  be  immediately  transformed  t  the 
other. 

3.  A  CLASS  OF  TOO  PARAMETER  FAMILIES  OF  DECREASING 
FAILURE  RATE  DISTRIBUTIONS 

We  shall  postulate  a  class  0  of  concave  hazard  functions  which  has 
special  properties  encompassing  both  mixed  exponential  and  "work  hardened” 
life  distributions. 

Each  element  of  this  class  will  generate  a  two  parameter  family,  which 
is  a  subset  of  the  DFR  distributions.  Reliability  will  be  of  the  form, 
where  Q  €  ^  is  a  given  hazard  function, 

t. 

R(t)  =  e-^^  t  >,0;  a, (3  >  0, 

Here  q  =Q'  is  a  decreasing  hazard  rate  which  is  twice  differentiable, 

standardized,  i.e.,  q(o)  =  1,  and  it  and  the  induced  functions  \|c  and  £,  where 

<p(x)  ■  xq(x)  and  c(x)  =  1  +  xq'(x)/q(x)  (3.1) 

satisfy  _  , 

,o  . 

1  ip  is  increasing 
2°  q  is  log-convex 

3°  £  has  a  limit  at  ®  and  0  _<  s(®)  _<  1  f 
alternatively,  sometimes  the  stronger  condition  is  assumed,  viz., 

3'  C  is  decreasing. 

The  question  arises,  'lihere  did  such  assumptions  come  from  and  what  dis¬ 
tributions,  if  any,  satisfy  them?”  One  sees  immediately  that  if  q  decreases  then 

v 

q  decreases  for  any  y  €  (0,1)  so  that  D  is  closed  under  fractional  powers  C 
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of  the  hazard  rate.  The  Afanasiev  generalization  of  the  lorax  distribution 
has  hazard  rate  qA  given  by  (2.3).  One  sees  that  qA  is  decreasing 
and  lnqA  is  convex.  Moreover,  one  checks  easily  that  4^  is  increasing  and 
in  this  case 

SA(t)  =  [1  +  (l-r)t]/(l+t)  .  (3.2) 

is  decreasing  with  £A  («)  =  1-y.  So  in  this  archtypical  example  assumptions 
1°,  2°,  and  3'  are  met. 


In  the  Bernoulli  mixture  of  the  two  exponentials  we  find  from  (2.5) 


-  5  - 


rv 


(e-vX+r) 


2  —  vx 

,  q£(x)  = 


(e‘vX+r)2 


(3.3) 


and  since  r  >  1, 


q£(x)  =  v^re  vX(r-e~vX)/(e  vX+rP  >  0, 


thus 


_ _ 2  -vx  ' 

Sg(x)  =  1 - 


(Xe“vx+r)  (e-vx+r) 


(3.4) 


so  we  see  that  £g(«)  =  1  and  that  3°  is  satisfied. 

We  next  show  that  tg  is  an  increasing  function.  To  see  this  we  note 
that  tg(x)  -  0  for  all  x  >  0  if,  after  simplification, 

Xe~vx  +  r2evx  2  v2rp.x  -  r(l  +  X) ,  for  all  x  >  0  . 

Let  y  =  evx  then  the  inequality  above  becomes  clearly  true  for  1  <y  syQ, 
where  in  yQ  5  (1  +  X)/(  X-  1).  Let  y=yQt,  for  t  >  1,  then  the  inequality  to 


be  proved  becomes 
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- —s —  +  r2y„  2  vrp-  ,  for  all  t  >  1. 

.2  o  t 

But  the  right  hand  side  is  maximized  at  t  =  e,  so  it  is  sufficient  to  have 
yQ  >  vu/re,  which  is  implied  by 


P  s 


X-l 


1  +  (X-l)e2/(X-1)  X(X  +  1) 

(The  approximation  is  always  less  than  the  bound.) 


(3.5) 


We  new  must  show  that  2rq  is  convex.  To  prove  this  it  is  sufficient  ; 
to  show  that  QgQg  -  (lg)  .  But  substitution  and  simplification  from  (3-3) 
and  (3-4)  into  the  equation  above  shows  (we  spare  the  reader  these  details) 
that  this  inequality  is  true  for  all  x  >  0  iff  r  >  ,  if  and  only  if 

p  <  (1  +  A)'1  (3.6) 

We  claim  that  (3-5)  and  (3.6)  would  be  virtually  always  true  in  practi  i 
since  if,  for'  example,  X  should  be  as  high  as  10,  then  the  probability  of 
passing  such  a  bad  component,  with  a  failure  rate  ten  times  the  nominal 
design  rate,  must  not  exceed  .082.  This  would  seem  to  be  a  reasonable 
assumption,  at  least  for  firms  that  intend  to  remain  in  business.  Furthermot 
this  demonstration  shows  that  not  all  mixed  exponential  (and  hence  not  all 
DFR)  distributions  satisfy  our  three  assumptions  and  that  there  are  DFR 
distributions,  which  are  not  mixed  exponential,  which  do . 


Let  us  consider  the  failure  rate  q  of  the  gairnia  distribution  itself,  <■- 
assuming  a  shape  parameter  0  <  v  <  1  and  unit  scale  parameter.  For  x  >  0, 
the  failure  rate  satisfies  the  relation 

03  oo  M 

l/q(x)  =  f  (1  +  ~)v  ^e-tdt  =  x  f  (1  +y)v-^e_y  dy  (3.7)j 

q)  "0 

- - - _ - - -  - -  - ~ - ”■' - 1. 


Prom  the  first  equality  we  see  1/q  is  increasing  for  any  0  <  v  <  1.  Prom 
the  second  we  see  [xq(x)]  ^  is  decreasing  for  this  DPR  distribution;  thus 
1°  is  satisfied.  Again  using  the  first  equality  in  (3.7)  we  find  that 
£(x)  is  non-negative,  does  not  exceed  unity  and  approaches  one  as  x  -*■  ~  . 
Hence  3°  is  satisfied. 

To  check  that  2°  is  satisfied  also  is  more  difficult;  we  examine 

00 

£nq(x)  =  (v-l)ln  x  -  x  -Zn[  tv_1e“tdt]  . 

It  is  thus  sufficient  to  show  that 


Unq)’’(x)  =  — I— v  +  q’(x)  >  0  for  all  x  >  0. 


(3.8) 


From  (3.7)  we  have,  after  differentiation, 


/  (1  +  i)v-2 1  e-tdt. 


q2(x) 


Thus  to  prove  (3-5)  it  is  sufficient  to  show  that  l-vixq'(x), 
which  is  the  same  as 

Kvr  •-**]’ *  fW1 

Making  the  substitution  y  =  t/x,  then  writing  the  squared  term  as  the  product 
of  two  integrals,  one  in  u  and  the  other  in  v,  the  inequality  becomes 


oo 

i tl  + 


vV-1  -ux, 
u)  e  du 


/ (1  .  v)"'1  e'^dv  >£  (1  .  y)v'2  ye'^dy. 


Making  the  change  of  variable  u  +  v  =  y  on  the  left,  the  inequality 
becomes 

oo  y  oo 

f  /  (l-u^+y+  uy)v  *du  e  yxdy  >  J  (1  +  y)v  ^  ye  ^xdy. 


0  0 
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Thus,  it  is  sufficient  to  show  that  the  integral  on  the  left  exceeds  that 


on  the  right,  i.e.. 


Jr  2-|V-l 

J  1  +  ^jr  l  y  du  >  T  for  all  0  <  u  <  y  <  «, 

o  L  “  j  y 


which  is  clearly  true  .  [~I 


Of  course  not  all  DFR  distributions  can  be  obtained  by  an  exponential 
mixture  such  as  given  in  (2.0).  For  an  exact  description  of  the  extreme 
points  of  the  class  of  DFR  distributions,  see  Langberg  et  al  [9].  Moreover, 
only  particular  mixing  distributions,  thought  to  be  of  practical  interest, 
will  concern  us  here,  along  with  their  corresponding  Afanas'ev  generalizations. 
Thus  we  will  examine  only  a  subclass  of  the  DFR  distributions. 

In  a  recent  paper,  McNolty,  Doyle  and  Hansen  (1980)  have  dealt  with 
the  mixing  problem  by  examining  some  mathematical  methods  for  inverting  a 
general  mixed-exponential  reliability  to  obtain  the  mixing  density.  They 
have  also  discussed  some  of  the  physical  interpretations  of  the  relationship 
involved.  We  attack  a  related  problem. 

We  now  ask  what  are  sufficient  conditions  on  the  mixing  distribution  so 
that  assumptions  1°,  2°  and  3°  will  be  satisfied  for  any  Afanas'ev  general¬ 
izations  of  mixed -exponential?  In  the  following  discussion  we  always  omit 
the  limits  of  integration  when  they  extend  from  0  to  «. 

Theorem  1:  If  the  mixing  distribution  G  is  such  that  the  function 
K(x,y)  =  G(y/x)  for  x,y  >  0 

is  totally  positive  of  order  2,  then  1°  is  satisfied,  i.e.,  ^ 

is  increasing. 


( 
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Proof:  Making  a  change  of  variable  in  the  definition  we  see 

t(t)  =  J ye~ydG(y/t)  //•  ydG(y/t) . 

Upon  integrating  numerator  and  denominator  by  parts,  we  obtain 

j G(y/t)ye“ydy 

t(t)  +  1  =  — -  • 

J G(y/t)e-ydy 

Let  t^  >  t2,  we  must  show  \|e(t^)  -  'yhis  true 

jG(y/t1)ye'ydy,  jG( y/t-L)e  ydy 

'  °  S  /  r 

jG(y/t2)ye'ydy,  J  G(y/t2)e  ydy 


By  applying  the  basic  composition  formula  of  Karlin,  see  e.g.,  p.  100, 
Barlow  and  Proschan  (1975) ,  to  the  right-hand  side  above,  it  becomes  equal  to 

r  r  GfyfV’  G(W  yie‘n •  e'71 

II  X  ^1^2- 

J  J  G(y1/t2),  G(y2/t2)  y2e  Y2,  e  72 

yl  <  y2 

Clearly  the  second  determinant  is  negative  and  that  the  first  is  negative, 
by  definition  of  TP-2,  can  be  seen  by  setting  =  t^,  =  t^.  [] 

The  hypothesis  of  this  theorem  has  a  relation  with  Polya-Frequency 
functions  of  order  2,  if  it  were  expressed  as  the  difference  of  the  logarithms 
rather  than  as  a  ratio. 
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•  1 


As  we  have  previously  seen  for  a  binomial  mixture  of  exponentials,  some 
restrictions  were  necessary,  namely  p  <  (1  +  VX)-\  to  insure  that  toq  was 
convex  everywhere .  We  now  give  some  conditions  that  a  mixed  failure  rate 
be  log-convex  on  (0,  «). 

Theorem  2 :  If  for  a  distribution  G  on  (0,°°)  we  define  the  convolution 
(2) 

G  =G  *  G  and  the  related  function 

A  ^ 

(t-  u)c  dG(t)  dG(u)  for  >  0  (3-9) 

J 

u=0  t=0 

then  the  failure  rate  of  the  G-mixed  exponential  distribution  is 
log-convex  iff  for  every  exponential  random  variable  Z  we  have 


g2U) 


-r  r 


* 


where 


E  2AG(Z)  >  EZa'q  (Z), 

AQ(y)  =  J  4  G2(y-4)  d  G(2)(.i)  . 


0 

Proof:  From  equation  (2.1)  we  find  for  fixed  t  >  0  that 

2  2  2 
(^M-q)  Unq)"  =  ^  |a2M’0+|Jl'1 


where 


X1  e~Xtd  G(X)  for  i  =  0,1, < 


Thus  we  may  write 


(3.10) 


M-0M,3_Mll‘L2  =y x3g  XtdG^x^  /e  ytdG(y)  -  f  x2e-xtdG(x)  J  ye_ytdG(y) 


Setting,  for  notational  simplicity, 

d  H(x,y)  =  e”^x+y^tdG( x)dG(y)  , 

we  find  it  can  be  written 

.2 


// 


x  (x-y)  d  H(x,y)  +  JJ  x  (x-y)  d  H(x,y), 
x>y  x<y 


■1 
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and  by  making  an  interchange  of  variables  in  the  second  integral,  obtain 


■// 


M-qM-3  -  =  /  I  ( x+y ) ( x-y )  dH(x,y), 

x>y 


By  a  similar  argument  we  obtain 


•// 

x>y 


(x-yT  dH(x,y ) . 


From  arguments  of  syrmetry  we  see  that 

t^0  =  JJ  x  dH(x,y)  =  JJ  y  dH(x,y) 

and  so 

^0M'1+M'l  =  Jf  ^x2+xy^  dH(x,y)  =  JJ  (v2+xy)  dK(x,y) 


Hence  we  see  that  for  c  =  2 C ^  s  with  the  obvious  change  of  variable 

c(£nq)"  =  JJ  (u+v  1  dH(u.,v) 


x>y 


?  / 

( x+y ) C  x-y )  dK(x,y) 


(u+v)2  dH(u,v)  II  (x-y)2  H(x,y) . 


Mow  let  w  =  u+v  ,  and  simplify  to  find  that 


2c(£nq)  -  J  [  JJ  (x-y)2  (x+y-w)  dH(x,y)]  w  e  twdG^2\w). 
lantity  in  square  brackets  above  can  be  written 

[...]=  J  J  ( x-y ) 2  ( x+y-w)e-t  ^ x+y  ^  dG(y )  dG( x) . 


x=0  y=0 


Letting  x  +  y  =  4  we  obtain,  in  the  case  the  density  G'=g  exists, 

c. 


/oo  ✓•co 

J 


(4-2y)2(4-w)e-t4g(4-y)  g(y)dydi. 


0  y=0 

For  this  case  we  set 


-I 


0 


g2(4)  =  /  U-2y)  g(4-y)  g(y)  dy  for  4>0, 

and  verify  that  G2  =  g2>  In  the  general  case  the  quantity  becomes 

[...]=  J  U-w^’^dG^U) 


Hence  for  some  c^  >  0  we  have 


~  =  J  J  wU-w)  e_t(4+^  dG2(4)  dG(2)(w).  (3.11) 


c^(^nq) 


(2) 

Assuming,  only  for  notational  convenience,  that  both  g2  and  g  exist 
we  have 


c^( £nq)‘ 


•f 


-tx 


/: 


(2) 

( X— 4  )  ( 2  4-x)  g0U)  gV  (x-4;  C4 


Let  us  set 


f(4,x)  =  4  g2(x-4)  g^U) 


for  0  <  4  <  x, 


(3. 


then  we  can  rewrite  the  quantity  in  braces  above,  after  breaking  the "integral 

into  two  parts  and  changing  variables  in  the  second,  as 
x/2 

{...}=/  (x-24)[f(4,x)-f(x-4,x)]d4  =  2Aq(x)  -  XAq(x) . 

0 

The  second  equality  is  obtained  by  integrating  each  term  of  the  difference 
by  parts  and  simplifying  the  resulting  expressions  where 

AQ(y)  =  J  f‘  4g2(x-i)g^2^  (4)  didx. 

Jo  J  0 


By  utilizing  the  Liebniz  rule  for  change  of  the  order  of  integration 
this  can  be  seen  to  be  equal  to  the  expression  for  AQ  given  in  the  hy 
The  proof  of  the  theorem  is  completed  upon  noting  that 


c1(2nq)~  =  J  e"tx[2AG(x)-xA'  (x)]  dx  .  Q 


We  new  give  some  sufficient  conditions  that  £nq  is  convex. 


(3. 
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Corollary  1:  If  the  mixing  distribution  G  is  such  that  the  induced  function 
defined  in  (3.10)  is  concave  increasing  then 

2Aq(x)  2  xA£(x)  for  all  x  2  0 

and  the  failure  rate  is  log-convex. 

The  proof  is  obvious  since  is  positive,  decreasing. 


Corollary  2:  If  the  mixing  distribution  G  has  a  density  g  =  G'  and  the  ratio 

x  g^(x)/g2(x)  for  x  >  0 

(2) 

is  monotone  decreasing,  where  g2  and  gv  are  the  corresponding 
derivatives  of  the  functions  defined  in  theorem  2,  then  the  failure 
rate  is  log-convex. 

Proof:  It  is  sufficient  to  note  that 

f(4,x)  2  f ( x— a , x )  0  <  4  <  x/2, 

where  f  was  defined  in  (3.12).  Thus  we  must  show 

4g  (.6)  g2(x-*)  2  (x-4)g2(,j)  gv  '  0  <i><  r/2 , 

which  upon  division  is  seen  to  be  guaranteed  by  the  monotone  behavior  of  the 
ratio  assumed  in  the  hypothesis  .  jj 

Remark:  If  G  is  a  f(a,  1)  distribution  then  one  verifies  directly  that 


x  g(2)(x)  =  2a  +  1 

g2(x)  x 


for  x  >  0 


which  decreases. 


Ey  Corollary  2  we  are  assured  that  the  gamma  mixture  of  exponentials 
will  have  a  failure  rate  which  is  log-convex:  A  fact  which  can  be  verified 
directly  from  (1.4). 
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We  now  examine  the  behavior  of  the  induced  function  £  where  Q  e  0, 
We  note  that  always,  since  q'  <  0  and  >  0,  that 


1 


>  1  + 


xo/ 

q 


V_ 

q 


>  o 


If  X,  is  a  given  function  bounded  between  zero  and  one,  differentiable  at 

zero  with  £(Q)  =  1  and  £'(0)  2  0,  having  a  limit  at  »  and  for  which  the 

ratio  £-1  is  increasing  for  t  >  0  then  by  regarding  the  expression 
t 


(3.1.3)  as  a  differential  equation  in  the  unknown  function  q  the  solution 

q(t)  =  exp  {-J 


dx}  for  t  >0 


defines  a  decreasing  failure  rate,  which  is  log-convex,  standardized  and 
for  which  t  q(t)  is  increasing,  i.e.,  it  satisfies  the  assumptions 
1°,  2°  and  3°.  We  note  that  £(t)  =  1-y  +  y(7  +  t)-1  generates  the 
Afanas'ev  generalization  of  the  gamma-mixed  exponential  model.  We  now  study 
the  behavior  of  the  induce.1  h  in  the  case  of  a  general  mixed  exponential 
in  particular  its  asymptotic  behavior  as  determined  by  the  behavior  of 
the  appropriate  mixing  distribution  G  near  0. 


We  do  not  consider  as  reasonable  (for  our  application)  mixing 

distributions  which  are  of  infinite  order  at  the  origin,  such  as 

_x-2 

G(x)  =  e  for  x  >  0.  We  admit  for  consideration  only  distributions  which 
are  either 

(i)  discrete  in  some  neighborhood  of  zero  or 


(ii)  of  finite  order  at  the  origin,  viz.,  there  exists 
<  >0  such  that  x" K  G(x)  ■+  c>0  as  x  -*•  0 
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Theorem  3:  For  any  G*mixed  exponential  distribution  the  induced  function 
rr  has  the  properties  that  £r(0)  =  1,  is  initially  decreasing,  bounded 
above  by  unity,  and  the  limit  exists  with 

=  1,  if  G  is  of  type  (i)  and 
£q(“)  =  0,  if  G  is  of  type  (ii),  moreover 
>  0  iff  is  increasing. 

Proof:  Considering  equation  (2.1),  we  can  write  for  any  t  >0, 

,  -  “tG  (t)  _  t  /  e'Xtd  G2(x)  -  (3. 1*0 

1-^t)  ^tr~  r  .yt  (2) 

J  ye  ^rd  G^  '  (y) 

using  the  methods  and  notation  of  theorem  3  . 

Thus  we  see  from  (3-1.4)  that  £  (o)  =  1,  and  since  the  right  hand  side 
is  positive,  that,  1  2  £(t).  We  see  £  initially  decreases  linearly  with  slope 
q'(o)  since 

^to)  =  lim  £(x)-l  _  ,,  \ 

x-Vo  x  4  v  /- 

Making  use  of  classical  Tauberian  theorems  on  LaPlace  transforms, 
e.g.,  Widder  (1946)  p.  181,  we  obtain  the  limiting  behavior  of  at  ». 

The  last  claim  follows  from  the  identity  Kq  =  'I'q/Qq  •  D 

We  now  consider  conditions  on  the  mixing  distribution  which  will  insure 
that  1  -  K  is  either  monotone  increasing  or  unimodal.  Such  behavior  can 
often  be  easily  checked  in  each  specific  instance  but  we  have  a  sufficient 


condition  in 
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Theorem  4:  If  the  associated  function  Aq,  as  determined  from  the  mixing 
distribution  G  in  (3.10),  is  such  that  the  kernel 

K(x,y)  =  AG(y/x)  for  x,y  >  0 
is  totally  positive  of  order  2,  then  the  function  1  - 

U 

induced  by  the  failure  rate  of  the  G-mixed  exponential 
distribution  will  have  at  most  one  mode  in  [0,  <*>] 

Proof:  By  taking  the  derivative  of  equation  (3.14)  and  considering  only 
numerator  we  find,  utilizing  notation  from  the  hypothesis  of  theorem  2,  that 

num[-S'(t)]  =  JJ  C-ty(x-y)  +  y]  e~(x+y)t  dG(2)(y)  dG2(x). 


By  comparing  the  first  term  above  with  equation  (3.11)  and  using  the 
representation  in  (3-13)j  3lS  well  as  making  a  change  of  variable  in  the  second 
term  to  obtain  a  convolution,  we  find  after  simplification 


•/ 


(y-2)e  y 


AG(y/t)  dy. 


We  note  the  function  (y-2)e  y  for  y  0  changes  sign  exactly  once, 

and  so  by  applying  the  variation  diminishing  properties  of  Polya  frequency 

functions,  see  e.g.,  p.  93,  Barlow  and  Proschan,  we  conclude  that  (t) 

G 

changes  sign  at  most  once.  Therefore,  1  -  possesses  at  most  one  mode.  [] 


In  any  life  length  model  one  is  interested  in  the  distribution  resulting 
when  independent  components,  from  the  same  family,  are  connected  in  a  series 
system.  Prom  the  concavity  of  the  hazard  Q  e  2  folwws  the 
Remark:  If  components  with  independent  life  lengths  (<^,8^  for  1  =  1,. 

for  some  Q  e  2,  are  in  series  then  the  life  of  the  system  T  will  satisfy  the 
stochastic  inequality 


3 


n> 


T  =  min  T,  >  Jn  (la.,  la.&./Ia.) 

i  ~  Q  1  11  1 


with  equality  when  p  =  p. 
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Another  property  one  would  ask  of  any  model  in  which  the  failure  rate  is 
initially  decreasing  is  what  improvement  can  be  made  by  "burning-in"  a  component 

having  such  a  life?  In  practice  it  is  often  assumed  that  as  a  result  of  a  bum-in 
period,  surviving  components  are  exponentially  lived.  In  fact,  bum-in  tests  are 
often  required  in  electronic  component  procurement  with  a  statement  of  the 
ultimate  failure  rate  so  obtained. 

Of  course,  not  all  decreasing  failure  rate  distributions  do  become  constant 
after  some  finite  initial  period,  but  that  is  an  assumption  which  is  often 
thought  to  be  appropriate.  This  indicates  the  importance  of  the  second  model 
introduced  and  of  its  utility  in  a  determination  of  the  economic  value  of  the 
stochastically  extended  life. 

The  residual  life  T^  of  a  component  with  new  life  length  T  and  a  bum-in 
of  duration  x  >  0  is  the  conditional  life  remaining  after  time  t  given  that 
it  is  alive  then;  i.e.,  T  =[T — r | T  >  x]. 

The  residual  life  of  any  G-mixed  exponential  is  again  a  mixed  exponential 
but  with  a  different  mixing  distribution.  The  residual  life  will  have  density 
f(y)  =  ce~^dG(Y),  where  c  is  the  normalizing  constant. 

But  is  the  class  H  closed  under  burn-in?  A  bum-in  of  duration  t  for  a 
component  with  hazard  Q  e  2  will  yield  a  residual  life  T^  with  hazard  rate 

q  (t)  =  q  (t  +  x)  for  t  >0. 

Clearly  q^  is  not  standardized  but  one  sees,  after  brief  reflection,  that  \|f 
is  increasing,  q^  is  log-convex,  £  approaches  a  limit  in  [0,1],  if  X,  does, 
and  lastly,  £  is  decreasing  If  X  Is.  Thus  the  answer  is  affirmative,  except 
for  standardization. 

It  is  easily  seen  that  Tx  is  stochastically  larger  than  T  for  all  x  >  0 
but  usually  the  burn-in  time  x  is  increased  until  any  incremental  decrease  in  the 
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residual  failure  rate  is  not  worth  the  incremental  cost  of  testing.  This  point 
will  necessarily  be  somewhat  different  for  each  particular  Q. 

Remark:  A  bum-in  of  t  units  of  time  on  a  component  with  a  new  life  T  ~  JA(cs,0) 
will  yield  a  residual  life 

Tt  ~  J^a(1  +  0/(1  +  t3)]. 

If  T  “  J  (a,g:p,\)  then  a  bum-in  length  x  will  only  alter  the  proportion  of  high 

D  ~ 

failure  rates,  viz., 

Tt  ~  JB(a,p:p',X) 

where  the  altered  proportion  is  given 

p-  =  pe-XT|!  /(pe-^S  +(l-p)e-T?). 


4.  ESTIMATION  OF  PARAMETERS  USING  INCOMPLETE  SAMPLES 

The  samples,  that  are  obtained  when  components  having  a  DFR  life 
distribution  are  tested,  are  virtually  always  incomplete  in  the  sense  that 
testing  is  stopped  before  all  components  have  failed.  A  datum  on  a  component 
that  "failure  has  not  yet  occurred  after  a  specified  life"  is  called  in  practice 
an  alive  time  or  a  run-out .  Samples  containing  such  observations  are  said  to 
be  cen4o-ted.  Samples  in  which  life  tests  are  truncated  at  some  preassigned  ordered 
observation  occur  infrequently,  in  our  experience,  when  electronic  components  are 
tested. 


It  is  assumed  throughout  this  section  that  we  are  given  a  sample  vector 
t  =  (t1,  ...,  t^,  ...,  t^)  where  t^,  ...,  t^  a re  ordered  observations  of  times 
k+1,  .  .  . ,  t  are  the  ordered  observed  alive-times,  with 


of  failure  while  t, 


1  <  k  <  n.  All  observations  are  presumed  to  have  been  obtained  by  testing 
components  having  a  J^(a,(3)  distribution  with  unknown  parameters  a  and  B,  but 
with  Q  e  2  and  Y  specified. 

We  now  introduce  notation  for  the  two  empirical  distributions  (call  them 
and  F  )  of  the  times  of  failure  and  of  the  sample,  respectively.  We  set 

Fk(y)  =  ^  (#  of  t.  <  y  for  i  =  1,  .  .  .  ,  k} 


and  similarly  for  F  ,  and  we  make  the  notational  convention,  to  be  used 
subsequently  for  any  function  g,  that  its  transform,  according  to  the  empirical 
distribution,  is  denoted  by  the  proper  affix, 

/do  ^ 

9(xt)dF  (t),  g ( x )  =  /  g(xt)dF,  ( t) . 

0  J  0 

Some  results  will  now  be  given  on  maximum  likelihood  estimation  of  the 
unknown  shape  and  scale  parameters  in  the  case  of  censored  samples  with 
Q  e  2  and  y  e  (0,1]  specified. 

The  sample  vector 

t  =  (t^,  .  .  tk,  .  .  .  tn)  for  1  <  k  <  n 
corresponds  to  the  observed  events 


[T^  =  t^]  for  i  *  1,  ....  k  and  [T^  >  t  ^  ]  for  i  =  k  +  1,  .  .  . ,  n. 

By  definition  the  log-likelihood,  after  substituting  from  (2.1),  is  given  by 

k  n 


k  «n(aB)  + 


y  £n  q(t18) 

i=l 


a 


Z  «(“ie  >  ■ 
'<=1 


Dividing  by  the  constant  k,  we  write  resulting  function  of  the  parameters 
a, 0,  given  the  vector  t_,as 
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L(a,  p|t)  =  Jin  a  +  2n  B  +  lnq(8)  -  Q(0),  (4.1) 

* 

where  we  have  made  use  of  the  rotational  convention  introduced  earlier. 

Parenthetically,  for  given  B  >  0  the  likelihood  L(*js,  t)  is  concave 
on  (0,  »)  and  the  MLE  of  a  exists  uniquely  and  can  always  be  obtained  from 
the  equation  L'(a|B.t_)  =  0. 

We  thus  have 

Remark:  When  the  scale  parameter  B  >  0  is  given,  there  exists  a  unique 
MLE  of  a.  say  a,  given  explicitly  by 

a  =  k/n  Q  (8).  •  (4.2) 

This  result  is  well  known.  If  8  is  known  and  Q  given,  concave  or  not, 
then  the  values  y^^  =  QCt^B)  for  i=l,  .  .  n  can  be  calculated.  They  are  the 
alive  and  dead  times  from  an  exponential  distribution  with  unknown  failure  rate 
a.  The  total  life  statistic  divided  by  the  number  of  failure  yields  the  usual 
maximum  likelihood  estimate  of  the  mean  life. 

It  is  also  true  that  not  any  set  of  n  positive  numbers  (t^,  .  .  t^) 

with  1  2  k  5  n  designated  as  failure  times  and  the  remainder  as  alive  times  can  be 
used. to  estimate  uniquely  both  the  unknown  parameters  for  any  Q.e  D.  In  some  sense 
the  sample  must  be  close  to  what  would  be  likely  from  such  a  hazard  function. 

We  examine  the  case  for  unknown  R  in 

Theorem  5 :  When  the  shape  parameter  a  is  presumed  known  and  a  sample  vector 

t  =  (tj,  .  .  t^,  •  •  • »  tR)  and  Q  e  Dare  given,  there  exists  an  MLE 

of  B,  denoted  by  8,  and  defined  implicitly  as  the  smallest  positive  root 
of  the  equation 

K  (x)  -  t  (x)  =  0 


(4.3) 
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only  when 

Inf  £t-1(y)  <  -2-  <1  (4.4) 

y  >  o 

Remark:  The  condition  (4.4)  is  determined  principally  by  the  behavior  of  Z, 
since  t  is  always  increasing  by  assumption  2°,  mapping  (0,  »)  into  (0,“  ). 

For  example, £  may  be  monotone  decreasing,  as  is  with  S^(")  =  1-Y,  while 
£  only  decreases  initially  but  is  concave-increasing,  ultimately,  with 

zB(~)  -  i- 

Proof  of  theorem  6:  From  (4.1)  we  see  the  likelihood  L(fl|a  ,  t)  can  be  written, 
neglecting  constants  not  depending  upon  B,  as 

L(fl|a,  t)  =  -  ~y-  Q(B) 

If  we  define  A(0)  =  L'(s|a,  t),  we  see  we  must  determine  the  roots  of 

A(x)  =  *(x)  -  -f-  7  (x) 

which  confirms  equation  (4.3)  as  the  appropriate  one.  To  check  the  necessary 
condition  we  see  the  likelihood  increases  initially  as  x  increases  since 
A(Q)  =  £(0)  =  1.  But  t  is  always  increasing,  while  Z  is  only  initially 
decreasing,  since 

lim  £'(0)  =  £'(0)  t  <  0. 
x+0 

But,  by  assumption  3°,  the  limit  Z  (»)  is  in  the  unit  interval.  Thus  for  a 
solution  to  the  equation 

£  C^rty)]  *  f'ory>  0 

to  exist,  the  range  of  the  composed  function  must  contain  the  value  an/k,  as 
stated  in  the  hypothesis.  [] 
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1 

| 

i 

\ 

We  see  that  the  second  smallest  positive  root  of  equation  (4.2),  if  it 

•  1 

exists,  will  be  a  minimum  likelihood  estimate  of  3.  This  situation  can  occur 
frequently.  If  there  are  either  more  than  one  local  maxima  to  L(p|a,  t_) ,  i.e.,  j 

three  or  more  roots  to  (4.2),  or  there  is  no  maximum  in  (0,«),  then  this  would  j 

indicate  that  the  presupposed  choice  of  a  and/or  of  Q  e  D  should  be  reexamined.  j 

That  is  to  say,  either  the  initial  failure  rate  a  is  not  of  the  right  magnitude 
to  reflect  the  nuirber  of  first  failures  observed  or  the  induced  function  £  does 
not  decrease  monotonely  over  a  sufficiently  long  interval. 


A  simpler  situation  exists  in  the  following  case. 

Theorem  6:  If  Q  e  2  is  an  Afanas'ev  generaliztion  of  a  G-mixed  exponential, 
when  G  is  of  type  (ii)  and  0  <  v  <  1,  with  a  known  then  for  a  given  sample 
the  NILE  of  3  exists  iff 


Proof:  We  note  by  theorem  3,  since  G  is  of  type  (ii),  that 
the  induced  function  for  the  G-mixed  exponential,  is  [0,1]. 
we  find  £  -  1-Y  +  so  that  the  range  of  %  is  [1-v,  1], 

U 

equation  (2.1)  that 


f(x)  =  x 


1-v 


/ G(y/x)d(ye  y) 
~~j  G(y/x)d(e_y) 


the  range  of  Zn, 
Letting  q  =  (qG)V 
We  find  from 


so  that  lim  t(x)  =  -  .  thus  by  equation  (4.4)  the  result  follows.  [] 

X— 


Corollary  3:  Under  the  hypothesis  of  theorem  7,  if  Kq  is  monotone  then  so  is  x. 
and  the  MLL  of  p  exists  uniquely,  i.e,,  there  is  at  most  only  one  solution 
to  the  equation  (4.3). 


t 
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We  now  turn  to  the  estimation  problem  when  both  a  and  6  are  unknown. 

Theorem  7:  For  a  given  sample  t,  with  Q  e  D  specified  and  a,B  both  unknown,  a 
MT.F  of  say  £,  exists  as  the  smallest  positive  root  of  the  equation 


* 


* 


£(x)  -  4>(x)  =  0  for  x  >  0 

Where  4>  =  \jT/Q  ,  iff  the  sample  satisfies  the  inequality 

2tt<7  c^.5) 

When  $  has  been  determined,  the  NILE  of  a,  say  is  then  given  by  an 
analogue  of  equation  (4.2),  namely, 

a  =  k/n  Q( g ) .  <^-6) 


Proof  of  theorem  7:  Consider  the  likelihood  function  L(a,3|t),  as  given  in 
(4 .1),  defined  over  the  positive  quadrant.  All  the  stationary  points,  which  are 
determined  by  t_,  can  be  found  from  the  simultaneous  solution  of  |^-  =  0,  =  0. 

This  yields  the  two  equations  in  a  and  8  ; 


K  (0)  =  +(Q);  4“  =  IT  Q(0)- 

* 

*  Combining  these  into  a  single  equation  in  the  unknown  Q,  we  are  led  to  seek  3  as 

a  zero  of  the  difference  £(x)  -  $(x)  for  x  >  0  where  4>  =  i/Q  .  If  we  maximize  the 
likelihood  with  respect  to  a  for  any  value  of  8  by  substituting  a  =  k/nC(G)  into 

*  (ti.l)  and  obtain  L(k/nQ(9),3|  t ) ,  then  neglect  constants  independent  of  3,  we  obtain 
a  function,  call  it  B,  which  can  be  written,  new  using  argument  x,  instead  of  3,  as 


I 
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B ( x )  =  2.nq(x)  -  An[Q(x)/x] . 

Interchanging  the  order  of  integration  in  Q(x)/x  we  obtain 


Q(x)/x  =  t  f  q(tx)  dW  ( t)  =  t  •  q*( x) 

J  0  n 


where  W  is  the  distribution,  with  density  given  by 
n 

W|;( t)  =  [1  -  Fn{t)]/t  for  t  >  0. 

Again  neglecting  constants,  we  have 


<•«,  oo 

B(x)  =  2.nq( tx)  dF.  (t)  -  UnM  q(tx)  dW  (t)J 

JO  '  Jo 


=  £nq(x)  -  £nq*(x)  . 


(4.7) 


Using  a  Maclaurin  expansion  of  q  we  see 


B(x)  =  x  q 1  (0)  [t  -  t*]  +  O(x^) 


as  x  •*  0, 


where 


-I 


-'-l 


t  =|  t  dF^(t)  and  t  =  {  t  dWn(t)  =  t  /2t  . 


Because  q'(0)  <  0,  it  follows  that  B  is  positive  in  a  neighborhood  of  zero 
iff  t  <  t*.  Since  B(0)  =  B(")  =  0,  we  must  have  B'(.x)  =  0  for  some  x  >  0. 

Moreover,  one  verifies  that  xE'  =  ^  -  4>  and  so  the  MLE  can  be  found  as  the 
smallest  positive  root  of  this  difference.  (The  second  smallest  root,  if  it  exists, 

A 

will  be  a  minimum  likelihood  estirate.)  Having  determined  B  one  uses  the  partial 

A 

derivative  equated  to  zero  to  obtain  a.  [] 

We  new  discuss  the  situation  when  the  sarple  fails  to  satisfy  the  condition 
2t  t  <  t2,  and  MIX'S  do  not  exist  uniquely.  Inis  means  that  the  model,  i.e.,  the 

choice  of  Q,  may  not  be  appropriate  and  either  a  constant  failure  rate  model  or  a 
convex  failure  rate  model  (one  that  is  initially  decreasing  and  ultimately  increasing) 
may  be  indicated  rather  tliar.  a  DF?  model  -with  3  e  o. 
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5.  THE  COMPUTATION  OF  G  FOR  CENSORED  SAMPLES 


As  a  matter  of  practical  calculation  we  are  concerned  with  the  smallest 
root  of  the  difference  £  -  4>,  in  the  case  when  £  is  monotone  decreasing. 
Equivalently,  let  us  consider  the  composite  function  f(x)  =%  [<t>(x)]  in  a 

neighborhood  of  zero  with  the  location  of  the  smallest  crossing,  if  there  is 
more  than  one,  of  the  ^5°  line. 


An  alternative  expression  for  $  is  $  =  xQ/Q,  with  Q  a  convex  function 

decreasing  between  t  and  t  •  q(«),  while  the  smoothed  Q  ,  viz,,  O/x,  decreases 

L  * 

between  the  same  limits  at  a  slower  rate.  It  follows  that  Q/x  >  Q  >  0  so  that 
0  <  $  <  1.  Then  $  begins  at  unity,  initially  decreases  at  a  decelerating  rate 
and  tends  ultimately  to  £(«) .  To  see  this  note 


f 


1  im  /Q 

CO 


Tim  |V-V(«t)N(t«k(tx)dFn(t) 
J  t1'V(xt)Vq(tx)dFn(t) 


since  by  Widder  loc  cit  there  exists  a  u  e['(\l]  such  that  xuq(x)  -*•  a  t  0,  as  x  -*  ». 
The  composite  function  f  =  x,  $  behaves  in  a  neighborhood  of  zero  as  a 

A 

contractive  map,  being  initially  greater  than  x,  then  crossing  at  B  and  then 
being  below  x  for  a  range,  perhaps,  thereafter. 

Thus  we  know  that  successive  iterates 

B1+i  =  f(Bi)  for  i  =  0,1,2  .  .  . 

will  converge  to  B  >  0  as  long  as  0q  <  B°,  the  next  larger  zero  of 
K  -<t>,  if  one  exists.  Otherwise  the  iteration  will  converge  to  zero. 


Moreover  for  the  special  cases  in  (3.2)  the  inverse  X 


can  be  easily  found. 


u 
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*  ■ 


i 


k'. 


Our  computational  procedure  based  on  the  sample  is  as  follows: 

Algorithm:  Given  t. ,  .  .  . ,  t,  as  failure  times  and  t.  t  as 

-  I5  *  k  k+1  5  n 

censored  live  times  from  a  DFR  distribution,  proceed  as  follows: 

~  —  ~~2 

(i)  Compute  the  sample  moments  t,  t,  t 

~~2  _ 

(ii)  If  t  <  2t  •  t,  assume  the  observations  are  from  an  exponential 
distribution  with  failure  rate  X  and  estimate  it  by 


nt 

~~2  _ 

(iii)  If  t  >  2t  t,  assume  the  observations  are  from  a  DFR 

distribution,  with  prescribed  Q  z  D  . 

Using  the  sample  functions  explicitly  given  by 

n  n  _  k 

$(x)  =  x  I  tiq(xti)/  E  Q(xt.)  C(x)  =  1  Z[tiq'(xt.)/q(xti)], 

we  guess  3Q  ,  then  iterate  using  the  inductive  step; 

given  B.. ,  compute  <t>(B.j)  and  calculate  such 

tha  t 

we  find  §  =  1 i m  ,  and  compute 

a  =  k/  Z  Q(t.B). 

1  1 

The  nature  of  the  intersection  guarantees  that  within  the  region  when 
$  and  5  both  decrease  the  iteration  will  rapidly  converge,  with  a  rea¬ 
sonable  first  guess.  When  the  functions  q  and  Q  are  simple  a  small  pro¬ 
grammable  electronic  calculator,  such  as  the  HP-67,  can  be  used  to  obtain 
these  estimates. 


f 
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We  now  present  some  data  sets  from  two  different  lots  of  flight  control 
electronic  packages.  Each  package  has  recorded,  in  minutes,  either  a  failure 
time  or  an  alive  time.  .An  alive  time  is  the  time  the  life  test  was  terminated 
with  the  package  still  functioning  and  is  denoted  by  an  affix  +. 

First  Data  Set 

X,  8,  10,  59+,  72+,  76+ ,  113+,  117+,  124+,  145+,  149+,  153+,  182+,  320+. 
Second  Data  Set 

37,  53,  60+,  64+,  66+,  70+,  72+,  96+,  123+. 

One  checks  that  both  data  sets  satisfy  condition  (4.5)  so  that  both 
parameters  can  be  estimated  in  a  gamma  mixed- exponential  model.  Then  using 
the  estimation  techniques  derived  previously  in  this  paper  we  have  the 
following  estimates: 

Data  Set  1  Data  Set  2 

a  =  . 0453, 0  =  1.03  a  =  .420,8  =  .01 

A  statistical  test' to  determine  whether  the  data  require  a  constant  or 
decreasing  failure  rate  was  run  on  the  data  from  sets  1  and  2.  For  data  set 
1  we  reject  constant  failure  rate  (in  favor  of  decreasing  failure  rate)  at 
the  .10  level.  For  data  set  we  cannot  reject  the  constant  failure  rate 
assumption.  In  this  case,  however,  the  constant  failure  rate  estimates  for 
reliability  and  the  mixed  exponential  estimates  for  reliability  are  close. 


7.  Conclusion 


We  are  primarily  concerned  with  the  problem  of  estimating  the  hazard 
rate  of  a  component,  such  as  an  integrated  circuit,  assuming  that  it  has 
either  a  constant  hazard  rate  or  a  decreasing  hazard  rate  of  specified 
functional  form  with  shape  and  scale  parameters  unknown.  The  influence  of 
the  data  is  different  in  this  case  than  in  the  more  usual  case  of  hazard 
rate  known  to  be  increasing. 

If  a  component  has  a  life  distribution  with  an  increasing  failure  rate, 
the  information  necessary  to  estimate  its  parameters  must  contain  failure 
times.  In  practice  this  means  that  if  there  are  few  observed  failures,  within 
a  fleet  of  orperational  components,  there  is  little  information  with  which  to 
assess  their  reliability.  If  a  component  has  a  constant  failure  rate  then 
both  failure  times  and  alive  times  contribute  equally  to  the  estimation  of  its 
reliability.  This  study  suggests  that  if  a  component  has  a  life  distribution 
with  decreasing  failure  rate  it  is  the  alive  times  within  the  data  which  con¬ 
tribute  principally  to  the  estimation  of  the  parameters  (and  thereby  to  the 
determination  of  the  reliability)  since  only  one  failure  observation  is  required 
even  to  estimate  two  parameters,  presuming  the  data  is  ample. 

Note  that  for  a  sample  of  size  two,  both  of  which  are  failure  observations, 

the  inequality  (4.5)  cannot  be  satisfied  since  t,^  >  0  implies  that 
—  2  2  2 

2(t)  >  (t^  +  )/2;  but  if  t^  is  a  failure  while  t^  is  an  alive  time, 

for  which  (1  +  /2)t-^  <  t^,  then  the  inequality  is  true. 

To  illustrate  the  typical  behavior  of  the  likelihood  B,  as  given  in 

equation  (4.7),  let  us  assume  that  the  sample  distributions  F,  and  W  as 

k  n 

defined  previously  satisfy  Fk  >_  W  ;  this  implies  t  >_  t  .  (Strictly  speaking 
this  condition  cannot  ever  be  met  since  even  if  the  failure  times  are 
stochastically  smaller  than  the  sample  containing  failure  and  alive  times 


together,  i.e.  F,  >  F  ,  still  we  would  have  0  =  F,  <  W  on  (0,t,  )•  But 
°  k  —  n  k  n  r 

~  it 

since  F^  <_  on  (t^,®)  we  would  have  t  <_  t  if  and  only  if 


00 

/ 


[Wn(t)  '  Fk(t)]  dt  1  V  [1  "  1/2n  t] 


Because  t^  is  very  small  in  practice  this  is  virtually  always  true.) 
From  this  assumption  would  follow  q  >_  q*  and  since  in  is  concave,  that 
in  q*  >  (Jin  q)*.  Hence  from  (4.7)  one  sees  that 


£n  q  -  In  q  <  B  <  to  q  ■  (£n  q)*. 

Now  one  notices  that  the  lower  bound  is  alwasy  negative  while  the  upper  bound 
is  always  positive.  But  near  the  origin  B  is  approximately  equal  to  the  upper 
bound  while  for  large  values  of  its  argument  B  is  approximately  equal  to  the 
lower  bound  and  for  intermediate  values  B  makes  a  transition  between. 

Thus  for  any  sample  in  which  the  failure  times  are  stochastically 

_  2 

smaller  than  the  combined  sample  times  and  for  which  2t  t  <  t  we  find  a 

maximum  likelihood  estimate  of  8  as  the  smallest  positive  root  and  a  minimum 

likelihood  estimate  of  8  as  the  second  smallest  root.  In  the  case 
_ 

2t  t  >  t  ,  which  could  occur  with  complete  failure  data,  we  would  have  a 
minimum  likelihood  estimate  as  the  smallest  positive  root  and  a  maximum 

A 

likelihood  estimate  at  S  =  0.  The  frequency  with  which  this  occurred  under 
various  mixtures  was  studied  by  Sunjata  (1974).  Maxima  occurring  at  such 
boundaries  were  also  observed  by  Davis  and  Feldstein  (1979)  in  their  study.  Of 
course  by  arbitrarily  grouping  failures  several  local  extrema  of  B  can  be  con¬ 
structed.  This  is  regarded  as  being  of  little  practical  significance. 

The  usual  justification  for  using  maximum  likelihood  estimates  is  due  to 
their  asymptotically  optimal  properties,  and  to  their  asymptotic  normality. 

The  problem  of  obtaining  the  usual  sampling  distributions  of  the  MLE's  of  the 
parameters  obtained  for  these  DFR  models  seems  to  be  difficult,  not  only 
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because  the  estimates  are  only  implicitly  defined,  but  also  because  samples 
are  virtually  always  censored.  Furthermore,  the  usual  proofs  for  the 
asymptotic  optimality  of  the  MLE's  may  not  apply  when  censoring  is  of  a 
general  type  and  when  only  sparse  failure  data  are  available.  A  useful 
asymptotic  theory  must  be  developed  for  censored  data  sets  of  which  the 
life  of  electronic  pacages  of  integrated  circuits  are  an  illustration. 
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